# Objective: Determine likely % of houses with protection prior to 1990.
# Table output: Pre/Post 1990 | N | Fire-ready (FR) roof | FR siding | FR eaves | FR vent-screens
# 2 rows of data

pacman::p_load(tidyverse, sf, cowplot, lubridate, kableExtra)

rm(list = ls())

source("code/globals.R")
dins <- st_read(file.path(INPUT_PUBLIC, "DamageData/DINS_2013_2020/All_Incidents_2013_2020_DINS_GDB.gdb/"))

dins <- dins %>% st_drop_geometry()

# Only look at single family homes
dins <- dins %>% filter(STRUCTURETYPE %in% c("Single Family Residence Multi Story", "Single Family Residence Single Story"))

# Only look at incidents with a good sample of non-destroyed homes
dins <- dins %>%
  mutate(undamaged = DAMAGE == "No Damage") %>%
  group_by(INCIDENTNUM) %>%
  mutate(
    total = n(),
    shareundamaged = sum(undamaged) / n()
  )

plot(year(dins$INCIDENTSTARTDATE), dins$shareundamaged) # more recent years have much better sampling of undamaged homes.

dins <- dins %>% filter(shareundamaged >= 0.4)


# Add indicator for built before 1990
dins <- dins %>% mutate(before1990 = YEARBUILT < 1990)

table(dins$ROOFCONSTRUCTION, dins$before1990)
table(dins$EXTERIORSIDING, dins$before1990)
table(dins$EAVES, dins$before1990)
table(dins$VENTSCREEN, dins$before1990)
table(dins$WINDOWPANE, dins$before1990)


# Determine fire-readiness based on roof, siding, eaves, vent screens
# In general, count "Other" and "Unknown" as NA
dins <- dins %>%
  mutate(
    fr_roof = case_when(
      grepl("(Wood)", ROOFCONSTRUCTION) ~ FALSE,
      grepl("(Concrete|Metal|Tile|Fire Resistant)", ROOFCONSTRUCTION) ~ TRUE
    ),
    fr_siding = case_when(
      grepl("(Wood|Combustible)", EXTERIORSIDING) ~ FALSE,
      grepl("(Metal|Brick|Vinyl|Fire Resistant|Ignition Resistant)", EXTERIORSIDING) ~ TRUE
    ),
    fr_eaves = case_when(
      grepl("(Enclosed|No Eaves)", EAVES) ~ TRUE,
      grepl("(Unenclosed)", EAVES) ~ FALSE
    ),
    fr_ventscreens = case_when(
      grepl("(Mesh Screen <= 1/8|No Vents)", VENTSCREEN) ~ TRUE,
      grepl("(Mesh Screen > 1/8|Unscreened)", VENTSCREEN) ~ FALSE
    ),
    fr_windows = case_when(
      grepl("(Multi Pane|No Windows)", WINDOWPANE) ~ TRUE,
      grepl("(Single Pane)", WINDOWPANE) ~ FALSE
    )
  )

# Make into long
lng <- dins %>%
  select(APN, YEARBUILT, before1990, fr_roof, fr_siding, fr_eaves, fr_ventscreens, fr_windows) %>%
  pivot_longer(
    cols = starts_with("fr_"),
    names_to = "feature"
  )

# Set feature order
lng <- lng %>%
  mutate(feature = factor(feature,
    levels = c("fr_roof", "fr_siding", "fr_eaves", "fr_ventscreens", "fr_windows"),
    labels = c("Roof", "Siding", "Eaves", "Vent screens", "Windows")
  ))

# Construct table
tbl <- lng %>%
  filter(feature != "Roof") %>% # Roof data are not information -- no mapping from materials to fire resistance rating (e.g., `asphalt' could be good or bad)
  filter(before1990 == TRUE) %>%
  arrange(feature) %>%
  group_by(feature) %>%
  summarize(
    N = n(),
    fr_true = sum(value == TRUE, na.rm = T),
    fr_false = sum(value == FALSE, na.rm = T),
    fr_unknown = sum(is.na(value))
  ) %>%
  mutate(
    p_observed = (fr_true / (fr_true + fr_false)) * 100,
    p_ready_lower = (fr_true / N) * 100,
    p_ready_upper = ((fr_true + fr_unknown) / N) * 100
  )

# Save as kbl

kbl(tbl %>% select(feature:p_observed),
  col.names = c("Feature", "Homes", "Yes", "No", "Unknown", "% Fire-Ready"),
  format = "latex",
  linesep = "",
  digits = 1,
  format.args = list(big.mark = ","),
  booktabs = T
) %>%
  add_header_above(header = c(" " = 2, "Fire-ready?" = 3, " " = 1)) %>%
  write("results/fire-ready.tex")

### CHECKS ###

# Other time periods

tbl2008 <- lng %>%
  filter(YEARBUILT >= 2008) %>%
  arrange(feature) %>%
  group_by(feature) %>%
  summarize(
    N = n(),
    fr_true = sum(value == TRUE, na.rm = T),
    fr_false = sum(value == FALSE, na.rm = T),
    fr_unknown = sum(is.na(value))
  ) %>%
  mutate(
    p_observed = (fr_true / (fr_true + fr_false)) * 100,
    p_ready_lower = (fr_true / N) * 100,
    p_ready_upper = ((fr_true + fr_unknown) / N) * 100
  )

tbl98_07 <- lng %>%
  filter(YEARBUILT >= 1998 & YEARBUILT <= 2007) %>%
  arrange(feature) %>%
  group_by(feature) %>%
  summarize(
    N = n(),
    fr_true = sum(value == TRUE, na.rm = T),
    fr_false = sum(value == FALSE, na.rm = T),
    fr_unknown = sum(is.na(value))
  ) %>%
  mutate(
    p_observed = (fr_true / (fr_true + fr_false)) * 100,
    p_ready_lower = (fr_true / N) * 100,
    p_ready_upper = ((fr_true + fr_unknown) / N) * 100
  )


# chars are much more less likely to be reported for undamaged homes
table(dins$ROOFCONSTRUCTION, dins$undamaged)
table(dins$EXTERIORSIDING, dins$undamaged)
table(dins$EAVES, dins$undamaged)
table(dins$VENTSCREEN, dins$undamaged)
table(dins$WINDOWPANE, dins$undamaged)
